Synchronization of Huy gens' clocks and the Poincare method 



We study two models of connected pendulum clocks synchronizing their oscillations, a phenomenon 
originally observed by Huygens. The oscillation angles are assumed to be small so that the pendulums 
arc modeled by harmonic oscillators, clock escapements arc modeled by the van dcr Pol terms. The 
mass ratio of the pendulum bobs to their casings is taken as a small parameter. Analytic conditions for 
existence and stability of synchronization regimes, and analytic expressions for their stable amplitudes 
and period corrections are derived using the Poincare theorem on existence of periodic solutions in 
autonomous quasi-linear systems. The anti-phase regime always exists and is stable under variation 
of the system parameters. The in-phase regime may exist and be stable, exist and be unstable, or 
not exist at all depending on parameter values. As the damping in the frame connecting the clocks 
is increased the in-phase stable amplitude and period are decreasing until the regime first destabilizes 
and then disappears. The results are most complete for the traditional three degrees of freedom model, 
where the clock casings and the frame are consolidated into a single mass. 
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Abstract 



1 Introduction 



In 1665 Huygens reported that two pendulum clocks attached to a common frame tended, after a while, 
to run in unison with their pendulums always pointing in mutually opposing directions. He eventually 
attributed this "sympathy" to imperceptible motions of the frame that transfer motion from one clock 
to the other. With some modifications, Huygens' experiments were repeated by Ellicott, Euler and 
Poisson, among others, who observed similar behavior, see descriptions in JT\. Original attempts to give 
a quantitative description of this sympathy were based on linear theory. Korteweg in a paper of 1905 [1] 
showed that the motion of the linearized system decomposes into three normal modes, two of which involve 
sizable oscillations of the frame. He then argued heuristically that damping in the frame siphons energy 
away from these two modes, eventually leaving only the mode responsible for the anti-phase synchrony. 
While an improvement over Huygens' intuitive explanation, Korteweg's theory left all non-linear effects 
of clock escapements and coupling out of his quantitative analysis. Nonetheless, upgraded variations of 
modal analysis appear in some recent approaches to the problem [21 [3] . 

By the end of the 19th century coordinated behavior in systems of interacting self-excited oscillators 
was observed in many other acoustic, electric and mechanical systems, the phenomenon later came to 
be known as synchronization. About the same time, motivated by applications to celestial mechanics, 
Poincare adapted the method of small parameter to finding periodic solutions to non- linear systems. 
Somewhat later, the Poincare method was applied to synchronization problems by Blekhman and others, 
see [4J and references therein, providing a first non-linear theory of the effect. Since the method requires 
equations of motion to be analytic in phase variables Blekhman used the simplest analytic model of a 
self-excited oscillator, the van der Pol model, to describe clock escapements. His analysis predicted not 
only the previously observed anti-phase synchronization, but also the in-phase synchronization, when the 
pendulums always point in the same direction. According to Blekhman, both regimes always co-exist 
in the same system being triggered by different initial conditions, he also reported observing both in 
experiments. 

In the last decade interest in Huygens' sympathy was renewed by the work of Bennett et al. pj, who 
took great care to reproduce Huygens' experiments and developed a new non-linear model to describe 
them. Their work builds on Korteweg's modal analysis, but equations of motion include terms modeling 
clock escapements with impulsive kicks. Unlike Blekhman's their model does not predict a stable in-phase 
synchronization, but it does predict asymptotic vanishing of oscillations in one of the clocks, a phenomenon 
they call beating death. This is a direct consequence of having a threshold angle for escapements to engage, 
if at some point the amplitude of oscillations falls below it energy loss forces oscillations to damp out. The 
authors even expressed puzzlement over Blekhman's results since neither they nor Huygens ever observed 
the in-phase synchronization. However, Czolczynski et al. [5] later confirmed co-existence of both regimes 
experimentally using modern clocks with non-impulsive escapements. Numerical simulations in ^ i7i| also 
show co-existence for some parameter values. The common trait in all three cases is that the escapements 
involved do not have an engagement threshold. 

We should also mention that the in-phase synchronization was observed by Pantaleone [8j and others 
[IDi after replacing pendulums with metronomes. Moreover, for small damping in the frame the in- 
phase synchronization was the only one observed, in contrast to Huygens' setup. The reason is that 
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metronomes swing with amplitudes too large to be approximated by harmonic oscillators, changing the 
situation entirely. Existence of the the in-phase synchronization for large amplitudes was proved in 
using the van der Pol model for escapements. 

Our main purpose in this paper is to apply the Poincare method to the Huygens system and to 
compare predictions it makes to known observations. Our application of the Poincare method builds on 
the original work of Blekhman As already mentioned, one obvious limitation of the method is that 
all terms in the equations have to be analytic, so realistic escapements with impulsive kicks can not be 
accomodated. However, van der Pol approximation appears to be rather accurate as long as the pendulum 
amplitudes stay above the engagement threshold, see [8]. On the other hand, the method offers several 
important benefits. First, there is a natural small parameter in the system, the mass ratio of pendulum 
bobs to the frame, which incorporates heavy clock casings. Second, one can take into account not only 
escapements bit also the non-linear coupling term, which remains after pendulum equations are partially 
linearized for small angles. In [2j and [6j this term had to be dropped because both approaches rely on 
the system's evolution being linear between escapement kicks. However, its magnitude is comparable to 
the linear terms kept. Finally, the Poincare method identifies synchronization regimes rigorously without 
heuristic pre-simplifications required for modal analysis, and gives analytic expressions for their stability 
ranges and amplitudes. 

One of the method's key predictions, namely co-existence of the in-phase and the anti-phase regimes, 
appears to be at variance with some observations. We will show in Sections |3j [4] that this prediction is 
not due to the method's limitations, but arises from treating damping in the frame as a small parameter. 
In Section [3] we assume both coupling of the pendulums to the frame and damping in the frame to be 
small. To avoid neglecting the second derivative in the frame's position, we transform the system into a 
first order form in a way different from that in [Ij. However, our conclusions are the same in this setting: 
the in-phase and the anti-phase regimes co-exist for the same values of system parameters. 

On the other hand, once damping in the frame is treated as a regular parameter in Section [4] one sees 
that the two regimes are not created equal. While the anti-phase synchronization remains robustly stable 
under variation of its values, the in-phase regime is only stable when the damping is very small relative 
to coupling. As the damping gets larger, it destabilizes and then disappears altogether. Although the 
Poincare theorem, as any asymptotic result, does not prescribe what "small" is, simulations in Section [5] 
suggest that in Huygens' experiments values of damping were outside of the co-existence range. In other 
words, damping in the frame is a bifurcation parameter and experimental values for it can be both below 
or above the bifurcation values. To illustrate the point, numerical simulations in Section [5] are compared 
to experimental results in [21 [5]. 

As in |2J , additional clues are provided by a preliminary linear analysis in Section [2] on the system ob- 
tained by disregarding all non-linear terms. Although synchronization is often believed to be a non-linear 
phenomenon, it turns out not to be entirely so. Even the linear system has the sum of pendulum angles 
converge to zero, which is an anti-phase synchronization of sorts. Of course, due to linearity stable am- 
plitudes depend on initial conditions, it is their convergence to invariant values that is a non-linear effect. 
Also, in the linear model pendulums started with exactly equal phases have their oscillations asymptoti- 
cally damp out. Although van der Pol terms alter this asymptotic outcome, numerical simulations show 
that the system may initially follow the linear pattern. 
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In addition to the traditional three degrees of freedom model, where the clock casings and the con- 
necting frame are modeled by a single mass, we also consider a four degrees of freedom system, where 
each casing is represented by a separate mass (Section [g]). In the linear approximation this system reduces 
to just three degrees of freedom because the angle difference between the pendulums decouples from the 
rest of the variables. Similar system but with a piecewise-linear escapement term was considered in [6], 
where existence of a stable anti-phase synchronization is proved for some parameter range. 



2 Equations of motion and linear synchronization 

We start with a general setup of n identical pendulums ('clocks') of mass m and length I attached to a 
rigid frame of mass M, which is restricted to move horizontally. In its turn, the frame is attached to 
a fixed support by a weightless spring with stiffness k and damping c. Figure [T] Pendulum angles with 
the vertical are denoted by 6i and position of the frame is given by x. Horizontal restriction reflects the 
original setup of Huygens, where clocks were hanging off a wooden beam resting on the backs of two 
chairs. Experiments show that behavior of the system changes substantially if non-horizontal motion is 
allowed, and synchronization may not occur at all [5]. Equations of motion can be derived using the 




m m 



Figure 1: Rigid frame with pendulums attached to a fixed support. 

Lagrangian approach, see e.g. |9j, 

mt^Oi + mgl sm.{9i) + mix cos{9i) = Di, ^ = 1, ■■■,n 

(M + nm)x + cx + kx + ml Oi cos{9i) = ml Y^^=i ^1 sin(6'j) 

Here Di := D{6i, 9i) denote driving torques provided by the clock escapement mechanisms. The Poincare 
method, which is our method of choice, requires all terms to be analytic in their variables. Following 
Blekhman [1], we choose the simplest analytic term that admits a limit cycle, the van der Pol term 
D{9,6) := e(7^ — 6'^)6. Here e measures the strength of the escapement and 7 is the critical angle, at 
which the escapement switches from boosting to damping. This model does not always work at small 
amplitudes due to engagement threshold and impulsive kicks in some escapements, but on average gives 
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a good approximation when the pendulums are close to their limit cycles [8j. We also assume that 
oscillation angles remain sufficiently small so that the pendulums can be approximated by harmonic 
oscillators. Taking sin^j « 9i and cos^j ~ 1 we get a partially linearized system 

{mPOi + mgWi + mix = Di, i = 1, n 

(M + nm)x + cx + kx + ral Yl^=i = Y17=i ^i^i ■ 

- 1 and introducing a dimensionless frame position variable y = y we get a 

dimensionless system 

yei + 6, + y = Fi, i=l,...,n 



D' C kli TTli 

where Fi = — a = ; — , il^ = — and /3 = . Note that the only remain- 

T^gl {M + nm)y/g/l {M + nm)g M + nm 

ing non-linearities are on the right hand side of ([T]), namely the van der Pol terms and the non- linear 

n 

coupling term j3 OfOi. Before proceeding with analysis of this system it will be instructive to investigate 

1=1 

a linear version obtained by simply dropping the right hand sides in ([T]). The resulting linear system 

p'i + ^i + y = 0, i = l,...,n 

\y + ay + n''y + PY:i=iOi = ^ 

is used in [5] to attempt a simple explanation of Huygens' observations. Asymptotic behavior of this 
system can be analyzed using a Lyapunov function. It turns out that the sum of pendulum angles 
asymptotically vanishes. 

Theorem 1. Suppose that all coefficients in ^ are non-negative and 7 > 0. Let 9i = 9i{t) be solutions 
to it with arbitrary initial conditions, then 61 + ■ ■ ■ + 9^ > 0. 

i— >oo 

Proof. Set 9 := 9i + ■ ■ ■ + 9n, then adding the first n equations in ^ we get 

{9 + 9 + ny = Q 

\y + ay + n'^y + p9 = Q. 

Define the Lyapunov function by 

E ■= (3(9^ + 9^) + n(f]2y2 + ^ 2(3n9y = (39^ + (/3 - n(3^)9^ + nfl'^y^ + n{(39 + y)^. 

Then £^ > as long as (3 — > or < /? < ^. Given the definition of /3 above this holds even with 
strict inequalities. Moreover, if < /3 < - then ii^ = if and only if 9 = 9 = y = y = 0, i.e. E is positive 
definite. On the other hand, by ([s]) 

E := 2(39{9 + 9) + 2ny{y + fl'^y) + 2f3n{9y + 9y) = -2Pn9y - 2nif39y + 2-ff) + 2pn{9y + 9y) = -4njf. 



5 



By the LaSalle invariance principle any trajectory of (|3]) asymptotically converges to a trajectory contained 
in £""^(0) since E < 0. But for 7 > any latter trajectory must have y = 0, hence y = and y = const. 
The system then implies 9 = —6 = const and 9 = 6 = 0, so any trajectory within £^"^(0) does not leave 
the origin. Therefore, y > and 9 = 9i + ■ ■ ■ + 9n > 0. □ 

Now let us focus on the case of two pendulums that we are most interested in. 

Corollary. Let n = 2 in Theorem^ and set 5{t) := 9i{t) — 92{t). Then asymptotically 
(9,it) ~ sin(t + <P) ^ ^ I . 

Proof. Subtracting the equation for ^2 in Q from the one for 9i we have j + (5 = 0. Therefore, 5{t) = 
Asm.{t + (p) with A and (j) ^ in the statement of the corollary. Since asymptotically 9i{t) ~ —92{t) by 
Theorem[l|it follows that 9i{t) ~ y{t) 92{t). □ 



0.3 -I 




is measured in the number of cycles of a decoupled harmonic oscillator with the same parameters, (a) 
9i{t) and 92{t) (indistinguishable at this resolution) superimposed, (b) 9i{t) + 92{t). Linear pattern is 
followed for about 150 cycles. 



Thus, pendulum angles asymptotically align in anti-phase for any initial values. The difference with 
true synchronization is that the stable amplitude A depends on initial values. In particular, if the 
pendulums are started at exactly in-phase, 6{0) = 5(0) = 0, then A = and they will both asymptotically 
damp out. Even if initial phases are not exactly the same the asymptotic amplitude can be very small. 
Numerical simulations show that for large a and small coupling /? the non-linear system may follow this 
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linear pattern for many cycles before the van der Pol terms break the trend, see Figure [2j If an impulsive 
escapement is used instead, and the engagement threshold for it is above this asymptotic amplitude, one 
can expect that one or both pendulums will stop beating in finite time. 



3 Small damping in the frame 

Poincare method relies on non-linear terms entering equations being multiplied by a small parameter. In 
practice, this can be interpreted as some dimensionless parameter of the system being 'sufficiently' small. 

TTi 

In our case, it is natural to assume that /3 = — is small since masses of the pendulum bobs m are 

M + nm 

much smaller than the consolidated mass M of the clock casings and the frame. We will also assume that 

the escapement strength e := of the van der Pol term Fi = £(7^ — 0f)9i is small for the purposes of 

mgl 

the Poincare theorem (which may in fact be large in practical terms). Eliminating ij from the pendulum 
equations in ([T]) we get 



{9, + 9^-ay-n^y = F,-j^F, i = l,...,n; 

\y + ay + n^y = j^F, 

where 

n 

F := -n{ay + n^y) + (^^ + efOi - Fi 
1=1 

This transformation is needed to reduce the system to a first order form. In system Q 

a m 
P M+nm ^ 



(4) 



^ M+nm 



is a convenient choice for the small parameter. Since there can only be one small parameter we shall 
assume that e = fia and set Fi := 0(7^ — 6f)6i, so that Fi = fj.Fi. The system transforms into 



+ 9i - ay - n-^y = fi{Fi - F) , i = l,...,n; 



1 9 (5) 

[y + ay + n^y = fxF , 

and its Poincare analysis for n = 2 will be our main goal in the next section. 

In this section however, we wish to consider a simplified case when a is also considered small. We 

will often refer to o" = ; — as "damping in the frame", although it is also obviously affected 

{M + nm)^Jg/l 

by lenghts and masses of the pendulums and mass of the frame. Aside from serving as a stepping stone 
to the general case, the case of small a compares directly to the treatment in Chapter 5 of [4J. Similar 
extra assumption is made there, but system ([T| is not transformed into Q, making it necessary to neglect 
the backreaction of the frame on the pendulums given by y. In line with the above, we set a = fib and 
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subsume the y term in ([s]) into the right hand side. After neglecting the terms of order fi^ we finally 
arrive at 



1, 



(6) 



with Tj as above and 



I ^. + 0. _ f]2y ^ ^ ( _ j:) ^ 

n 
i=l 

System ^ can be rewritten in the first order form as x = Ax+/i$(x), where x = (^i, ^i, . . . , 6n,y, vY' ■ 
From now on we only consider the case of two pendulums, n = 2. Then explicitly 

1 
-1 
1 
0-10 
1 
-if 

The corresponding linear generating system (see Appendix for the terminology) is obtained by setting (j, = 
in ([g]). Poincare theorem gives sufficient conditions for existence and stability of periodic solutions to the 
original non- linear system that converge to periodic solutions of the generating system when — )■ 0. The 
latter can be found most easily by diagonalizing A = VAV^^. In our case, A = diag (i, i, —i, —i, —iQ,) 
and X = with ^ = [ifi, ip2, 11^1,1^2, u, v)'^ given by 



and 










T 









T 







J" 






1 n 



1 



2 



T {-iy + y) 



(7) 



u=\ {iVLy + y) , v = \ {-iVLy + y) 



Let us further assume that the system is non-resonant, i.e. the natural frequency of the frame Q is not an 
integer or a half- integer. In fact, it is realistic to expect \Q\ <C 1, meaning that the clocks oscillate much 
faster than the frame. Making the corresponding substitution in non-linear terms, namely 



we get the diagonalized system 



f(e) = v-ici>(vo, 



^ = Ae + /.f(o. 



The general solution to the generating system is ^ = A^ is ipk = Ofce**, tpk = /3fee 



u 



Ce'^\ V = De 



(8) 



(9) 



It is not periodic since Q may be irrational, and, in any case, the periodic solutions we are interested 
in should reduce to independent oscillations of the pendulums when the coupling is removed, i.e. when 
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/i = 0. To put it differently, we are looking for periodic solutions with the amplitude of frame motion 
of order at most fi. This means that we are looking for non-linear counterparts of the solutions with 
C = D = 0. 

Physically, parameter measures the strength of the non-linear coupling of the pendulums to the 
frame. The generating system describes two pendulums attached to an immovable frame, so they are 
decoupled. Once the frame is allowed to move the coupling is turned on and the behavior changes. It turns 
out that the pendulums eventually settle into oscillating synchronously with limit amplitudes determined 
by their escapement mechanisms. Depending on the initial phases this settled motion can be in-phase 
or anti-phase. In the anti-phase regime the frame remains almost motionless once the synchronization 
occured, the pendulums behave as if they were nearly decoupled. More sizable motions are required to 
sustain the anti-phase regime, speeding up the oscillations compared to decoupled pendulums. The next 
theorem makes this description more precise providing explicit formulas for limit amplitudes and periods 
to the first order in fi. It also implies that both regimes reemerge after small disturbances. 

Theorem 2. For small fi > 0, a,b > and 2il a non-integer system ([g]) with n = 2 admits both in-phase 
and anti-phase synchronized periodic solutions with limit amplitude 27, where j ^ is the critical angle 
in the van der Pol escapement. Aside from the trivial one, these are the only solutions that converge when 
IJ, ^ to pendulums independently oscillating with equal amplitudes, and motionless frame. The period of 
the anti-phase solution is T{n) = 27r-|- o(/i), and that of the in-phase solution is r(/i) = 2-k — |^^2 fJ' + o{n). 
Both solutions are stable. 



Proof. In this proof we use notation and terminology associated with the Poincare theorem, see Appendix. 
Recall that the diagonalized matrix of the generating system is A = diag {i, i, —i, —i, ifl, —iQ). We want 
our periodic solutions to reduce to oscillations of decoupled pendulums for ^ = meaning that only 



ai, . . . , 04 in ( A.l ) should be non-zero. Therefore, we must choose the first four eigenvalues to form the 
leading special group. The remaining two then belong to non-special groups, and the secondary special 
group is empty by assumption about $7. 



Our next task is to simplify constraint (A.3) from the Poincare theorem for our case. By definition, 
^ = {ifi, ip2, "01) ''P2, u, v)'^ and the relevant components of ^ are the first four. Since x = must be a real 
solution one can see from ([T]) that 4>k and ipk must be complex conjugates of each other for A; = 1, 2. In 
view of (A.l) we have ak+2 = ctk-, and we represent ak = r^e*'^''. Here is the half-amplitude of the kih. 
pendulum and (j)k is its phase. One can set (i)2 = ^ without loss of generality since time can be shifted by 
a constant, then synchronization is characterized by a single phase (pi = 4>. By assumption about equal 
amplitudes set ri = r2 = r. After some computations one can reduce (A.3) to a system of equations for 
r and (j), namely 







7^ (1 + 7^) sin(0) = 

^2 (1 + ^2) (^^2^<t> _ 1) = . 



(10) 



Discarding the trivial solution r = 0, the only remaining ones are r = 7, (/> = and r = 7, = vr, 
corresponding to the in-phase and the anti-phase synchronization respectively. Substituting the values of 
a. 



corresponding to each solution into (A.5) we get the first order period corrections. 
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Stability condition (A. 6) reduces to the following equations having only roots with negative real parts 



: (X + 07^) (^(1 - f)2)2^2 ^ ^^2^^ _ ^2-^2^ + (1 + 7^ 



0: 



vr : (x + 07^) ( (1 - n'^fx'^ + 07^(1 - 0^)^x + 87'' + 47^ + 1 



2\2. 



0. 



This condition is equivalent to all the coefficients of the above polynomials being positive, which happens 



if a > 0, 7 7^ and O 7^ ±1. Stability condition (A. 8) for non-special eigenvalues reduces to 6 > in both 

□ 



cases. 



Numerical simulations in Section [5] show that the limit amplitude usually gives a good approximation 
to actual amplitudes when is small. We believe that these are the only solutions to the Poincare system 



(A. 3) even without the assumption of equal amplitudes. However, we were unable to solve it for ri, r2 
and (f> without assuming that ri = r2. We can only show that if </> = or = tt is assumed in advance 
then ri = r2 follows from the equations. Based on simulations, we also believe that for small n these are 
the only attractors, i.e. the trivial solution is spurious (no periodic solutions with amplitudes of order /x), 
and no quasiperiodic or chaotic behavior occurs. 

We see that for small a there holds a counter-intuitive symmetry: both synchronization regimes 
always coexist and have the same stable amplitudes. This conclusion coincides with Blekhman's in [3] 
even though he neglects the second derivative of the frame's position and applies the Poincare theorem 
to system ([T]) directly. However, such coexistence of regimes was not observed in some experiments, see 
[2]. We shall see in the next section that this symmetry is a deceptive consequence of a being small, 
increasing it sufficiently eliminates the in-phase regime altogether. 



4 Damping in the frame and the in-phase synchronization 

In this section we no longer assume that damping in the frame is small and investigate how its magnitude 
influences synchronization. It is this assumption, it turns out, that puts the two synchronization regimes 
on an equal footing in Theorem [2] Recall that ^ was obtained from ([s]) after assuming that a is small. 
Without this assumption, after neglecting the terms of order /i^, we get instead 

ei + Oi-ay -n^y = fi{Ti-T) ^^^^ 
y + ay + Vi^y = 

with := a{j'^ — 6f)6i and 

n 

T := -n{ay + + J^(l + ^f)^. . 

1=1 
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As in the previous section, we rewrite this as a first order system x 
y, y)^. For two pendulums one has 



Ax + /i$(x), where x 






1 














-1 













a 











1 














-1 







a 

















1 
















—a 


is A 




diag 




-i, -i, 


1 

2 



and 










T 









T 







J" 





^ )j we omit 



the diagonahzing transformation V . 

Behavior of the system is most conveniently described in terms of a modified damping parameter 
o := —7- — -^T^ — For small values of and a it is close to the ratio of damping o to the escapement 

strength a. According to the theorem below the system undergoes two transitions as a and hence a 

2 

increase. For values smaller than a := 2(2+72) both synchronization regimes coexist and are stable. This 
agrees with the conclusions of the previous section. The settled frequency of the in-phase oscillations 
decreases as the damping in the frame grows, but still remains higher than for the anti-phase oscillations. 
Once the above value is exceeded the in-phase regime destabilizes, but remains in the picture until 

2 

the value o" := ^ is attained. At this point the in-phase regime disappears altogether. Parameter values 
in the experiments of Bennett et al. [2] (and presumably Huygens) seem to fall within this last range. 
This means that if clocks are started with nearly equal phases they will eventually settle into anti-phase 
oscillations with amplitudes determined by the escapement mechanism. This agrees with observations 
and indicates that a should not be treated as a small parameter when the Poincare method is applied. 



Theorem 3. For small fi > and a,a > system (11) with n = 2 admits anti-phase synchronized periodic 



solutions with the limit amplitude 27, where 7 7^ is the critical angle in the van der Pol escapement. 



In-phase synchronized periodic solutions can only exist if cr < \, where a :- 



'1-n' 



, and 



have amplitude 2\J \^2^ ■ ^^'^^^ from the trivial one, these are the only solutions that converge when 
— )• to pendulums independently oscillating with equal amplitudes, and motionless frame. In-phase 
solutions exist and are stable at least if the following inequalities hold: either aa < 1 and a < 



or 



aa > 1 and < a < 



7 



aa 2 ^ ^ 2(2+7^ ' period of the anti-phase solutions is T{fi) 
of the in-phase solutions is 



27r 



2(2+72)' 
o(/i), and that 



T{fi) = 27r 



(1 - J72)2 + 2 f + 



fl + 0(/i) 



Proof. The general outline of the proof is the same as in Theorem [2] The leading special group is again 
formed by the first four eigenvalues and the remaining two are non-critical for o" > 0. We wish to simplify 
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constraint (A. 3) from the Poincare theorem. After setting 

,(^2 = and (hi 



r^e*'^'' and taking into account that 



r2 



one can obtain: 



(7^ 

[i {n^ - 



r^) - 5e*<^ (1 + - a (l + r 



+ e- 



+2 



2) =0 

[i {n^-i)+a] 



'■2) a 



2a (l + r2) = 
= 0. 



Last equation immediately imphes that = or (/> = vr. Substituting these values into the first equation 
we get the desired half-amplitudes r for the in-phase and the anti-phase regimes. The second equation 
is then satisfied automatically. Expressions for periods are obtained by substituting the values of Ug 



corresponding to each solution into (A. 5) 



Stability condition (A.6) for = vr reduces to the equation 



(x + 07^) + a ((1 + 45)72 + 25) X + a-{l + 7^) (l + 72(0^ + 3) 







having only roots with negative real parts. For this to hold, all the coefficients must be positive. One can 
see by inspection that they are as long as a, a > and 7^0. The same stability condition for (p = is 
somewhat more involved: 



If < 



(x + 0(72 - 25)) ( x2 + 



2(2+7^) 



7^-25(2 + 7^) , 5(1 + 7^) , ^ ~ , 2. 



1 + 2aa 



-X + a 



£7(1 + 25)2 



1 + 2aaa + 7(1 — aa] 



0. 



and aa < 1, or aa > 1 and ""^ ^ ^ < o" < 



7 



2(2+72) 



then the coefficients in the second 



parentheses are positive. In both cases cr < ^ is automatically satisfied, so that the coefficient in the first 



parentheses is also positive. Other stability conditions ( A.7[A.8 ) do not apply here since the remaining 
eigenvalues are non-critical. □ 

Note that the no resonance condition on is no longer necessary because ct > suppresses any 
resonance. Moreover, setting = 0, i.e. eliminating the stiffness altogether, does not qualitatively alter 
the conclusions. In practice, when the escapement has an engagement threshold the in-phase regime 
disappears even sooner, at least when the asymptotic amplitude falls below the threshold. The in-phase 
regime is then precluded by one of the pendulums switching off completely. For this to happen, it is not 
even necessary for asymptotic amplitude to be below the threshold as long as the pendulum angles fall 
below it before the onset of synchronization. Numerical simulations (see Section [5]) also show that when 
the in-phase regime is unstable, a period of beats may precede the anti-phase synchronization. In the 
course of beating the two pendulums exchange much of their energy back and forth with pendulum angles 
falling far below their initial or eventual asymptotic values. This is likely the reason why vanishing of 
oscillations in one of the clocks was observed in the experiments of Bennett et al. |2J. 



12 



5 Numerical simulations 



In this section we use numerical simulations to illustrate the existence and stability conditions of Theorem 
[3j We also wish to demonstrate that the theory based on the Poincare method can explain, at least 




Figure 3: In-phase synchronization for parameters of Czolczynski et al. with ^i(O) = 0.25 and 02 (0) = 0.3. 
(a) 6*1 (t) (lighter) and 02{t) (darker) superimposed, (b) 9i{t) —92{t)- 

qualitatively, conflicting reports on the type of synchronization observed in experiments. Correspondingly, 
we tried to use the values of parameters from the experiments of Bennett et al. |2J and Czolczynski et al. 
[S]. Unfortunately, in the case of the former the information on these values is insufficient. 

First three pairs of graphs correspond to the values g = 9.81, m = 0.158, c = M = 11.856, k = 1.186, 
I = 0.269 from [5], we picked 7 = 0.122 for the critical angle and e = 5.047 for the escapement strength. 
Time is measured in the number of cycles of a decoupled harmonic oscillator with the same parameters. 
One easily checks that aa <\ and a < 2(2+^2) the notation from Theoremjsj hence both synchronization 
regimes exist and are locally stable. When the initial values are close (Figure |3]) the system synchronizes 
in-phase. As they are pulled apart, the asymptotic regime switches to anti-phase, but it is preceded by the 
depicted period of beats (Figure [4]). Finally, when the initial values are sufficiently far apart (Figure [5|, 
the anti-phase synchronization is achieved faster and without beats. This is in line with the observations 
of Czolczynski et al. In all three cases the stable amplitudes are nearly equal to the limit values given by 
Theorem [3l 

Next, we tried to simulate the experiments of [2j with g = 9.81, m = 0.082, / = 0.14. Unfortunately, 
there is no information on the values of stiffness k or damping c. The mass M was adjustable in ^ and 
we show results for two different choices of its value. We kept 7 = 0.122 and assigned the rest of the 
parameters as indicated in the captions. In the case of Figure |6] the second inequality from Theorem [3] 
fails and the in-phase regime is unstable. The system synchronizes in anti-phase but only after rather 
violent beating (individual oscillations are not visible on the figure because of a large time interval). The 
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Figure 5: Anti-phase synchronization for parameters of Czolczynski et al. with 9i{0) = 0.1 and ^2(0) = 
—0.3. (a) 9i{t) (hghter) and 02{t) (darker) superimposed, (b) 9i{t) + 92{t). 
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, with 



Figure 7: In-phase synchronization for parameters of Bennett et al., M = 6.143, c = k = 0.614, with 
di{0) = 0.1 and 6*2(0) = 0.3. (a) 6*1 (t) (lighter) and 6*2 (t) (darker) superimposed, (b) 6*1 (t) - 6*2(0. 
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onset of synchronization takes much longer than in the previous cases in hne with the observations of 
Bennett et al. In fact, we adjusted the missing parameters to match the onset times from [2j. Depending 
on the escapement threshold one may observe either the anti-phase synchronization or the vanishing of 
oscillations in one of the pendulums in this case. To stabilize the in-phase regime we had to decrease 
the damping by more than an order of magnitude (Figure [7|). We do not believe that such a value of c 
is realistic for the experiments of [2] (or of Huygens) explaining why they never observed the in-phase 
synchronization. One can see that the onset of the in-phase synchronization can also be preceded by 
beats. 

As one can see, Theorem [3] provides some theoretical guidance to observing coexistence of synchro- 
nization regimes, or lack thereof, in experimental setups. First, the van der Pol parameters a and 7 should 
be adjusted to match the system as closely as possible, e.g. along the lines of [8]. Then the inequalities 

2 

aa < 1 and a < 2(2+^^) give the parameter range, where both types of synchronization can be expected 

to manifest (the other pair of inequalities from Theorem ^ calls for unrealistically large values of damping 

2 2 — 

a) . When 2{2+'y^) ^ ^ ~2' expect beats to precede the onset of anti-phase synchrionization due 

to the presence of unstable in-phase regime. 



6 Frame model with two masses 

Starting with Korteweg most work on the Huygens' problem was based on variations of the three degrees 
of freedom model we considered above. In it the casings of the clocks and the connecting body (beam, 
platform) are consolidated into a single massive frame. However, it is the casings that make the frame 
massive, the connecting body itself is actually quite light [2l[5l|8]. In this section we consider a refinement 
of this model suggested by Dilao where casings are modeled by two separate masses M, and the 
connecting body is modeled by a massles spring of stiffness k and damping c, see FigjSj The rigid frame 
of the previous model (with mass 2M and detached from a fixed support) is recovered when k ^ 00. 
Each mass carries its own pendulum of length I and mass m, and, as before, the system is restricted to 
move horizontally. The positions of the casings are given by xi, X2 and the pendulum angles are 61, 92- 
Unlike Dilao, who uses a piecewise-linear model for escapements, we continue to use the van der Pol terms 
D{6,6) := e(7^ — 9'^)9 so that the Poincare method can be applied. Assuming small pendulum angles as 
before, the partially linearized equations of motion are 

ml'^9i + mgl 9i + mlxi = Di 

mP92 + mgl 92 + mlx2 = D2 

(M -|- m)xi + cxi — k {x2 — xi) + ml9i = ml 9\9i 

(M -|- m)x2 + CX2 + k {x2 — xi) + ml 92 = ml ^|^2 

with Di := D{9i,9i). Note, that Dilao drops the non-linear terms on the right in the last two equations. 
The reason is the same as for his choice of a piecewise-linear escapement. Dilao's approach requires 
evolution of the system to be linear as long as 9i stay below a critical angle, at which the escapment 
switches from boosting to damping. In contrast, the Poicare method only requires that the non- linear 
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Figure 8: Pendulums with clock casings connected by elastic frame. 



terms be analytic. Rescaling time as in Section [2j and introducing dimensionless parameters xi = lyi and 
X2 = ly2 we obtain a dimensionless system (compare to ([l])): 



Oi+ei + yi = Fi 
02 + 6*2 + 2/2 = 

yi + cryi - X (y2 - yi) + /3^'i = OlOi 
m + cry2 + >c {y2 - yi) + /3^2 = ^2^2 , 



where Fi 



Pi 

mgl 



mgl^l "il"^^ " - (Af+m)^' ^ ~ {M+m)g ""''^ " M+m ' 

parameter to be ;U := = ^ and set = as before. Separating terms multiplied by fi and 
discarding higher order terms we transform the system into the form 



kl 



and /3 



We choose our small 



61 + 01- ayi + K (y2 - yi) = A* f a(7^ 
^2 + ^2 - o-y2 - (y2 - yi) = /i f a(7^ 
yi + cryi - X (y2 - yi) = l^^i 

m + cry2 + X (y2 - yi) = /^-7^2 , 



9i) ^2 - -F2 



(12) 



where 



-^1 = -o-yi + >^ (2/2 - yi) + 6*1 ( 1 + and J='2=-cry2 - >c {y2 - yi) + 62 [I + Ol] . 
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System (12) can be rewritten in the first order form ^ = + where 
C = (6*1, 6*1, 6*2, ^2, yi, 2/1,2/2,2/2)^ and 
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+ Vcr^ - 8x ,-f7,0 



we 



omit the diagonalizing transformation. Note that one of the eigenvalues is zero, which indicates the 
presence of a rigid body mode. Indeed, one can see directly from Figj8]that the system as a whole is free 
to move rigidly. In the original model this was prevented by attaching the frame to a fixed support. But 
since this attachment, represented by fi'^ in (11), is a minor influence on synchronization we dropped it 
here for simplicity. As a compensation, we will only consider solutions that do not excite the rigid body 
mode. 

As the next theorem shows, qualitatively the behavior does not differ much from the case with the 
rigid frame. The in-phase regime eventually disappears, albeit at somewhat larger values of a. The effect 
of the frame stiffness x is similar in nature to that of the attachment stiffness O in Theorem [3l One 
should also keep in mind that M in the one mass frame model is twice its value in the two mass model 
here. This leads to an extra coefficient of 2 in amplitude formulas and ijl/2 in place of /i. With this in 
mind, the case of Theorem [3] for = is obtained in the limit x — )• 00, i.e. the limit of rigid frame. 



Theorem 4. For small n > and a, a > system (12) admits anti-phase synchronized periodic solutions 

(2x — 1)^ + fj^ ) , with amplitude 2^/ \^^°''^ o,nd period 



only if ^ > a an, where a, 



T{fj.) = 2-K 



(2x-l)(l + 72) ^ 
(2x- 1)V ? +a2 2 



+ o(m) • 



It admits in-phase synchronized periodic solutions only if 'j'^ > o'in, where o'in '■ — ^ 



2\-l 



with 



amplitude 2y \^^^" and period T[^) = 27r 
in the van der Pol escapement. 



1 + 7' 



1 + ^ + (j2 2 



+ o(/i) . As before, 7 7^ is the critical angle 



The proof is similar to that of Theorems [2] and [3] and we omit it. The only structural difference is that 
the eigenvalue has to be included into the special leading group, but the corresponding equation from 
(A. 3) is vacuously satisfied. Note the correspondence between ain and a from Theoremjsj Note also that 
Can < Cm for large x, meaning that the range of existence of the anti-phase regime is wider than of the 
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in-phase one, with disparity increasing as x — )• oo. Unfortunatley, we were unable to simplify the stability 
conditions to a tractable form in this case. But, by a continuity argument from Theorem [Sj at least for 
large x the in-phase regime should destabilize at smaller values of a than the anti-phase regime. 



7 Conclusions 

We used the Poincare method to study models of two connected clocks synchronizing their oscillations, 
a phenomenon originally observed by Huygens. The oscillation angles are assumed to be small so that 
the pendulums are modeled by harmonic oscillators, and the mass ratio fi of the pendulum bobs to their 
casings is taken as the small parameter of the Poincare method. Clock escapements are modeled by the 
van der Pol terms. The results are most complete for the three degrees of freedom model, where the 
connecting frame and the clock casings are consolidated into a single mass. The partial results on the 
four degrees of freedom model with casings represented by separate masses do not alter the qualitative 
picture. 

The anti-phase synchronization regime is dominant in the sense that it is universally stable under 
variation of the system parameters. In this regime (to the first order in n) the frame is motionless, and the 
stable amplitudes and periods of the pendulums are the same as for the decoupled van der Pol oscillators. 
The in-phase synchronization regime may exist and be stable, exist and be unstable, or not exist at 
all depending on parameter values. As the damping in the frame, or more precisely the dimensionless 



this regime first destabilizes and then disappears. Analytic conditions for existence and stability of 
synchronization regimes, and analytic expressions for their stable amplitudes and period corrections are 
derived based on the Poincare theorem. In particular, they provide explicit relations between escapement's 
strength and critical angle, and damping in the frame, that distinguish experimental setups supporting 
and not supporting the in-phase synchronization. 

Our results vindicate the intuitive explanation of Huygens' sympathy due to Korteweg, that the in- 
phase linear mode damps out faster than the anti-phase mode, because it requires more sizable motions 
of the frame, unless the damping is kept small on purpose. Interaction between the two modes explains 
the beats observed by Ellicott and others [1]. Beats are most likely to be observed when the in-phase 
regime is present but unstable. 

The following remarks are based on numerical simulations and are therefore somewhat speculative. 
When both regimes coexist, even if the in-phase one is unstable, the system may undergo a long period 
of beats before settling asymptotically. This is most likely to happen when the initial values are near the 
boundary between the basins of attraction of the regimes. In the course of beating the pendulum angles 
may fall significantly below their initial and asymptotic values. If escapement models with engagement 
threshold are used instead of the van der Pol terms the beats may cause one of the clocks to stop in 
finite time. When the damping is particularly heavy the system goes through a linear stage of anti-phase 
oscillations with asymptotic amplitude dependent on the initial values. If the engagement threshold is 
absent however, the clocks eventually synchronize in anti-phase with stable amplitudes. Thus, despite the 
limitations of the van der Pol model for escapements, we can account for all types of behavior observed 



parameter a 



c 



is increased the in-phase stable amplitude and period are decreasing until 
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in experiments [2| 15] and, to some extent, predict quantitatively when they occur. 



Appendix: Poincare method 

For convenience of the reader we briefly describe a method of studying periodic solutions to sytems of 
quasi-linear equations that we use in this paper. More detailed accounts can be found in [12], Ch.IX and 
[1], Ch. 10, proofs are given in [13j . Ch.V. 

Consider a first order autonomous quasi-linear system x = Ax-|-^$(x) and assume that the I x I matrix 
A can be diagonalized as A = VAV~^ by a nonsingular linear transformation V. Set F(,^) = V~^$(V^), 
then a change of variables x = Vy transforms the system into the form y = Ay -|- /iF(y). The linear 
system x = Ax, as well as y = Ay, is called the generating system. We are interested in the periodic 
solutions to the original system that converge to periodic solutions to the generating system when /i — )• 0. 
The latter exist only if some eigenvalues of A, the diagonal elements of A, are purely imaginary, and 
they are stable only if A has no eigenvalues with positive real parts. From now on we assume that 
all eigenvalues have non-positive real parts, and call purely imaginary ones critical. In synchronization 
problems one is typically interested in particular solutions to the generating system, which correspond to 
synchronous motion of decoupled oscillators. 

Let \ = iuj he a critical eigenvalue. Following Blekhman OH], we collect all eigenvalues of the form 
As = iusio with integer into the leading special group, and assume, without loss of generality, that they 
occur for s = 1, . . . ,k. Then the most general periodic solution of the generating system with period 
7" = 77 is given by 

0, s = k + l,...,l, 

where are complex amplitudes. According to the usual method of small parameter, one now looks for 
solutions y = y{t, /j,, a, P), with = y° + /3s for t = 0, to the original system as multivariable power 
series in fi,a,/3. If such solutions are periodic the period will be of the form T{fj,) = T{1 — 5{fj,)), where 
= 6iii + o(/i), since r(0) = T by assumption. 
Existence of periodic solutions to the original system now reduces to whether one can find power series 
that also satisfy 

y{T{fx),fi, a, p) - y(0, ^^, a, /?) = 0. (A.2) 

Indeed, by autonomy this implies that y is T(//) periodic. This is a system of / analytic equations with 
I + 1 unknowns /i, /3i, . . . , /?/, so one extra condition on /3's can be imposed. The system is solvable only 
if certain consistency constraints are satisfied. It turns out that ^ factors out of the first k equations 
and by setting their free terms to one obtains consistency constraints on 5i and the amplitudes ag. 



Eliminating 5i we finally end up with k — 1 equations (A. 3) for k amplitudes. Again, one extra condition 
can be imposed on a's, e.g. a^-i = a/c- This freedom reflects the possibility of shifting time in solutions 
to autonomous sytems. Thus, although the linear system has periodic solutions for arbitrary values of 
Os, non- linear periodic solutions only exist for some special collections of their values. Once the values of 



Oig are found from (A.3) the first order period correction 5i can also be determined, see (A.5). 
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If one is interested, as we are, in real-valued systems, then A is real and its eigenvalues come in 
complex conjugate pairs. The leading special group will contain both elements of a pair or none, so 
k = 2m and without loss of generality, Xm+s = for s = 1, . . . , m. Generating periodic solutions will 
also be real if and only if Um+s = ol- 

Stability of a periodic trajectory can be investigated by looking at the solution to the matrix equation 
Z = (A-|-/iV$(j;))Z with Z(0) = I, where V$(x) is the derivative of the non- linear term computed along 
the trajectory. This is the system of variations. Stability is determined by eigenvalues p{fi) of Z(T(/i)), 
e.g. it suffices that they were less than one by absolute value. For // = we obviously have Z{t) = e*"^, 
so Ps{0) = e''*^"^ and one can look for these eigenvalues in the form Ps{p) = e^''"^(l + xp + o{iJ,)). For 
non-critical eigenvalues |/0s(0)| < 1 and we automatically have |/0s(/^)| < 1 for small p by continuity. For 
critical eigenvalues we will have |/9<i(/i)| < 1 for small > as long as Re(x) < 0. 

For the purposes of stating the stability conditions it is convenient to distinguish further groups among 
the critical eigenvalues. Two eigenvalues belong to the same group if they are integer multiples of each 
other, or if they differ by inoj with n an integer. The point is that stability conditions decouple into 
separate conditions for each group. Recall that As = iHsUJ with integer Ug form the leading special group. 
The secondary special group is A^ = i{ns + ^)uj, and the remaining non-special groups have eigenvalues 

(r) 

of the form As = ivs , where r enumerates non-special groups. After heavy duty computations the 
conditions we described can be brought into a convenient form below, which is due to Blekhman [3]. 

Theorem (Poincare). Let y = Ay -|- ^F(y) be an autonomous quasi-linear system with diagonal A. 
Assume that all eigenvalues of A have non-positive real parts, and there is at least one purely imaginary 
eigenvalue X = iuj. For convenience, suppose that all eigenvalues of the form As = iUgUj with integer 
Ug are listed as first k. Then periodic solutions to the system becoming at p = periodic solutions 
to the generating system y = Ay with period T = ^ can only exist for such values of amplitudes 
ai, afc_2, Cik-i = ctkr that satisfy equations 

Qs{ai,..., Uk) = OkUkPs - OsUsPk = , s = l,...,fc-l, (A. 3) 

where 

Ps(«i,..., afc):=(Fse-'"^"*>:=;^ T Fs (aie*"^'^*, Q^e^"^"*, 0, O) e-™^"*dt , s = l,...,k. (A.4) 

^ Jo 

The first order period correction for these solutions is 



= (A^5) 



where a* are the values found from (A.3). 

These solutions will indeed exist and be stable (asymptotically orbitally stable) for small positive p if 
all the roots x of the following algebraic equations have negative real parts: 

- afenfc(5sjx = 0, s, j = 1, k - I ; (A.6) 
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{nj-ns)ojt \ _ ^ . 



JS] 



+ x 



0, s, J in the secondary special group; (A. 7) 



- Osj V XT 



- 0, s,j in the rth non-special group. (A. 8) 
// the real part of at least one root of these equations is positive then the corresponding solution is unstable. 
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